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SUMMARY 


As  part  of  a  combined  experimental  and  analytical  study 
of  vortex  shedding  at  sharp  edges  characteristic  of  strakes 
on  fighter  aircraft,  numerical  simulations  were  made  of  a  two- 
dimensional  and  a  three-dimensional  case.  The  former  is  a 
sharp-edge  flat  plate  normal  to  the  flow,  for  which  some 
experimental  data  were  taken  earlier  in  this  program.  The 
latter  is  a  double-delta,  sharp-edge  wing  at  angle  of  attack. 
In  both  cases,  a  free-stream  Mach  number  of  0.5  was  used.  The 
objective  of  the  work  was  to  investigate  the  suitability  of 
using  the  numerical  simulation  to  augment  experimental  flow 
field  data,  after  the  major  features  of  the  simulation  were 
verified  by  data. 

Two  existing  Navier-Stokes  codes  developed  at  Ames 
Research  Center  NASA  were  used  for  the  numerical  simulations. 
Neither  code  had  been  used  in  an  application  involving  flow 
separation,  sharp  edges,  and  the  aggregation  of  vorticity 
into  small  regions  in  the  flow  such  as  the  two  present  appli¬ 
cations  demonstrate.  The  applications  of  the  codes  were 
largely  successful,  in  that  for  both  applications  the  numeri¬ 
cal  flow  fields  generally  reproduced  the  gross  flow  features 
exhibited  by  experimental  data. 

The  larger  objective  of  demonstrating  the  suitability 
of  using  numerical  simulations  to  augment  experimental  data 
was  not  achieved.  The  principal  difficulty  in  both  applica¬ 
tions  was  mesh  resolution.  In  the  two-dimensional  case,  the 
mesh  resolution  was  inadequate  to  determine  the  details  of  the 
boundary  layer  on  the  windward  side,  particularly  near  the 
edge  where  experimental  data  is  very  difficult  to  obtain.  In 
the  three-dimensional  case,  the  problem  was  exacerbated  by  the 
requirement  to  store  more  information.  In  this  case  there 
were  differences  in  the  location  of  the  peak  suction  pressure 


on  the  wing  and  the  vorticity  distribution  above  the  wing 
which  were  caused,  at  least  in  part,  by  mesh  limitations.  The 
objective  is  felt  to  be  attractive  and  feasible,  but  consider¬ 
able  care  must  be  taken  in  the  numerical  simulation  to  permit 
the  mesh  to  capture  and  accurately  represent  important  flow 
phenomena . 
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1 .  INTRODUCTION 


The  use  of  strakes,  or  leading  edge  extensions,  on  recent 
fighter  type  aircraft  has  been  shown  to  improve  greatly  the 
high  angle  of  attack  aerodynamic  characteristics  of  these  air¬ 
craft.  The  beneficial  effects  stem  from  the  formation  of 
leading  edge  separation  vortices  on  the  highly  swept  strake 
and  the  subsequent  interaction  of  these  vortices  with  the  flow 
over  the  wing  to  produce  significant  lift  on  the  wing-strake 
combination  at  considerably  higher  angles  of  attack  than  is 
the  case  for  the  wing  alone.  These  benefits  have  been  derived 
primarily  by  cut  and  try  methods  in  wind  tunnel  tests  during 
aircraft  development  programs.  Some  studies  of  strake  aero¬ 
dynamics  have  been  undertaken  but  the  understanding  of  the 
fundamental  fluid  flow  phenomena  necessary  to  do  a  priori 
design  of  strake-wing-body  combination  does  not  exist.  The 
purpose  of  this  report  is  to  describe  a  fundamental  investiga¬ 
tion  undertaken  to  help  fill  this  need.  This  investigation 
builds  on  the  recent  work  that  has  been  started  on  strake  flows 
and  on  recent  advances  in  numerical  techniques  for  predicting 
fluid  flows. 

The  present  study  is  part  of  a  larger  program  conducted 
by  Nielsen  Engineering  &  Research,  Inc.  (NEAR)  for  the  Office 
of  Naval  Research  under  Contract  N00014-78-C-0388 .  The  first 
part  of  that  work  involved  some  two-dimensional  tests  in  a  wind 
tunnel  at  the  NASA  Ames  Research  Center  as  a  first  step  in 
developing  a  method  for  predicting  the  rate  at  which  vorticity 
is  shed  on  a  strake  edge.  That  work  is  reported  in  reference  1 
The  current  work  reported  herein  involves  numerical  simulation 
of  flows  over  surfaces  with  sharp  edges.  In  addition  to  this 
numerical  simulation  work,  the  program  was  to  include  some  wind 
tunnel  tests  on  a  highly  swept  wing  to  obtain  three-dimensional 


data  which  would  be  combined  with  the  numerical  and  two- 
dimensional  results  to  develop  and  evaluate  an  engineering 
prediction  method  for  leading  edge  shedding  rate. 

One  of  the  difficulties  in  using  the  three-dimensional 
data  is  the  fact  that  the  detailed  flow  field  data  acquisition 
is  a  time  consuming  process  and  the  quantity  of  data  that  can 
be  obtained  in  a  reasonable  tunnel  entry  is  not  as  great  as 
one  would  like.  Also,  it  is  very  difficult  to  obtain  boundary 
layer  measurements,  particularly  near  the  edge  where  the  flow 
separates,  which  is  the  region  of  greatest  interest.  There¬ 
fore  in  this  study  a  recently  developed  computational  tech¬ 
nique  for  solving  the  Navier-Stokes  equations  is  extended  to 
treat  flows  with  leading  edge  separation.  These  numerical 
results  can  then  be  used  together  with  the  experimental  data 
to  validate  the  theoretical  flows  and  to  define  the  flow  field 
in  greater  detail,  particularly  in  the  windward  side  boundary 
layer  and  near  the  edge.  The  combined  theoretical  (numerical) 
and  experimental  flow  fields  can  then  be  used  to  develop  the 
engineering  design  methods. 

Thus  the  driving  force  for  using  a  numerical  simulation 
of  the  flow  to  supplement  the  experimental  work  is  the  relative 
difficulty  of  obtaining  sufficient  detailed  flow  characteris¬ 
tics  from  the  latter.  The  kinds  of  details  required  for  the 
development  of  an  engineering  prediction  method  are  the  vor- 
ticity  in  the  subsequent  separated  vortex.  Some  of  the  exist¬ 
ing  Navier-Stokes  computer  codes,  especially  those  developed 
by  the  Computational  Fluid  Dynamics  (CFD)  Branch  at  NASA  Ames 
Research  Center,  are  available  to  help  define  the  flow  in 
more  detail.  However,  before  embarking  on  a  discussion  of  the 
numerical  compuations  it  is  helpful  to  be  more  specific  about 
the  problem  to  be  solved. 


Modern  strake  configurations  are  curved,  highly  swept  and 
blend  into  the  wing  to  give  a  "double  delta"  appearance.  To 
ensure  a  clean  separation  of  the  strake  vortex,  the  strake 
leading  edge  is  sharp.  The  free-stream  Mach  number  at  high 
angles  of  attack  where  strake  vortex  formation  is  important 
tends  to  be  subsonic.  Finally,  the  strake-wing  combination  will 
operate  at  high  angles  of  attack,  e.g.,  30°.  These  c  •'racter- 
istics  are  not  common  in  the  usual  computational  flu  dynamics 
problems.  The  only  available  three-dimensional  code  jlica- 
ble  to  the  problem  is  that  of  Lomax-Pulliam  of  the  C!  Branch 
at  Ames  Research  Center.  The  code  uses  a  fully  impl.  .  delta 
form  algorithm  developed  by  Beam  and  Warming  (ref.  2) .  The 
main  conceptual  difficulty  of  the  code  is  that  it  solves  a  "thin 
layer"  approximation  to  the  Navier-Stokes  equations.  This 
approximation  essentially  consists  of  the  retention  only  of 
those  second  derivatives  normal  to  the  surface.  This  allows 
treatment  of  separated  flows. 

The  remainder  of  this  report  is  concerned  with  the  two 
problems  studied.  The  first  problem  is  a  two-dimensional 
variant  of  the  configuration  examined  experimentally  in  com¬ 
pleted  work  for  ONR  (ref.  1) .  This  work  generated  information 
on  the  required  grid  clustering,  convergence  rates,  boundary 
conditions  at  the  sharp  edge,  and  the  ability  (and  limitations) 
of  the  code  to  predict  the  flow.  The  two-dimensional  experi¬ 
mental  results  of  reference  1  are  used  for  comparisons.  The 
information  gained  during  this  phase  of  the  investigation  was 
then  used  in  the  solution  of  the  more  complex  three-dimensional 
problem.  Since  the  two  problems  are  solved  by  different  codes 
and  are  sufficiently  dissimilar,  they  will  be  discussed  in  the 
two  following  sections  separately. 


2.  TWO-DIMENSIONAL  CASE 


The  problem  of  interest  is  the  two-dimensional  flow 
approaching  a  flat  plate,  as  shown  in  the  left  sketch  below. 

In  the  numerical  simulation,  the  plate  is  treated  as  infinitely 


thin  and  thus  has  sharp  edges.  For  comparison  purposes,  data 
are  available  from  reference  1  for  a  configuration  shown  in 
the  right  sketch.  The  plate  has  a  sharp  edge  with  a  finite 
angle  (a  cusp)  and  a  splitter  plate  extending  well  downstream 
to  prevent  alternate  vortex  shedding  which  can  occur  when  the 
two  sides  of  the  flow  can  communicate  with  one  another. 

The  data  (ref.  1)  which  are  available  for  comparison  are 
static  pressure  distributions  on  the  front  face,  two  orthogonal 
velocity  components  in  the  flow  field  near  the  edge  of  the 
plate,  and  the  location  of  the  reattachment  point  at  the  rear 
of  the  separation  bubble  behind  the  plate.  The  velocity  data 
can  be  used  to  calculate  vorticity,  rate  of  transport  of 
vorticity,  and  circulation  in  specified  contours  in  the  flow 
field. 


The  experiment  was  conducted  at  Mach  numbers  of  0.25  and 
0.5  at  a  Reynolds  number  based  on  the  half  face  width  in  the 
range  of  0.15  -  0.25x10®.  Thus  the  boundary  layer  on  the 
front  face  approaching  the  edge  is  laminar. 

2.1  Governing  Equations 

The  equations  governing  the  flow  adjacent  to  the  two- 
dimensional  strake  are  the  compressible  Navier-Stokes  equations 
Since  the  Reynolds  numbers  are  0.5x10®  or  less,  laminar  flows 
only  will  be  considered.  The  equations  are  in  Cartesian 
coordinates 
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Here 

Y  =  specific  heat  ratio 
p  =  pressure 
a  =  sound  speed 
Re  =  Reynolds  number 
<  =  thermal  conductivity  coefficient 
X  =  -(2/3)p,  by  Stokes  hypothesis 
p  =  viscosity  coefficient 
Pr  =  Prandtl  number 

The  velocity  components  are  nondimensionalized  with 

respect  to  the  free-stream  sound  speed  a^,  the  density  p 

with  respect  to  poo/  and  the  total  energy  e  with  respect  to 
2 

pooaoo  .  The  pressure  is  referenced  with  respect  to  yp,*- 

Since  the  two-dimensional  strake  is  a  thin  flat  plate, 
the  Cartesian  mesh  is  ideally  suited  for  this  case.  However, 
clustering  is  necessary  to  adequately  resolve  the  flow  field 
near  the  plate.  Thus  two  one-dimensional  clustering  functions 
are  introduced , 


S  =  Ux), 

n  =  n  (y) ,  (3) 

t  =  t. 


Under  this  transformation  the  governing  equations  become 
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The  transformation  Jacobian  is 


J  =  l/(x^yn)  (5) 

The  viscous  terms  are  resolved  only  in  a  thin  layer  near 
the  body  due  to  lack  of  sufficient  computer  capacity.  The 
viscous  terms  in  y  only  are  retained  and  those  along  £  (along 
the  strake  surface)  are  dropped  from  the  equations.  Thus 
equation  (4)  is  simplified  to 
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The  nomenclature  is  shown  in  Figure  1. 

The  thin  layer  model  is  a  useful  simplification  incorpor¬ 
ated  in  the  numerical  codes  used  for  this  study  (refs.  2  and  3) 


For  this  two-dimensional  case  where  £  is  normal  to  ri,  the 
viscous  terms  dropped  are  the  same  as  those  deleted  in  boundary 
layer  theory.  However ,  the  normal  momentum  equation  is 
retained  in  the  thin  layer  model  and  the  pressure  can  vary 
through  the  viscous  layer.  Therefore,  the  thin  layer  model 
does  not  incur  the  matching  problem  of  boundary  layer  theory 
with  the  inviscid  field  and  separation  points  are  not  singular 
points. 

On  the  body  surface  [n(y)  =  0  and  £(k)  <_  1],  the  no  slip 
condition  requires 


The  pressure  on  the  body  surface  is  obtained  from  the  normal 
momentum  equation  which  for  this  configuration  is  the 
y-momentum  equation. 

The  other  four  surfaces  of  the  flow  domain  of  Figure  1 
are  the  symmetry  plane,  inflow  and  outflow  planes,  and  the 
far  field  boundary.  Symmetry  conditions  are  imposed  on  the 
plane  of  symmetry.  Free~ stream  conditions  are  imposed  on  the 
inflow  plane  and  the  far  field  boundary. 

On  the  outflow  plane  zero  gradients  of  all  the  flow 
variables  are  imposed. 

2.2  Numerical  Algorithm 

The  numerical  algorithm  used  to  solve  the  transformed 
equations  and  boundary  conditions  is  the  time-implicit  numeri¬ 
cal  algorithm  developed  by  Beam  and  Warming  (ref.  2).  Since 
the  scheme  and  the  code  has  been  reported  elsewhere  (refs.  2 
and  3)  it  will  not  be  presented  here  again.  The  reader  is 
referred  to  the  cited  references. 


2.3  Checks  on  the  Numerical  Results 

The  numerical  results  were  obtained  on  a  thin  flat  plate 
of  unit  semispan  normal  to  a  flow  with  free-stream  Mach  number 
of  0.5.  The  Reynolds  number  based  on  free-stream  velocity 
and  semispan  is  0.5*10®.  The  schematic  of  the  flow  geometry 
is  shown  in  Figure  1.  The  inflow  and  outflow  boundaries  are 
at  eight  and  twenty  semi spans  upstream  and  downstream  of  the 
plate/  respectively.  The  centerline  is  a  plane  of  symmetry 
and  the  far  field  boundary  is  nine  semispans  from  the  center- 
line.  The  boundaries  of  this  domain  were  determined  by 
numerical  experimentation  so  that  the  flow  field  in  the 
neighborhood  of  the  plate  did  not  change  with  any  further 
extensions  of  the  boundaries. 

Checks  were  made  to  evaluate  accuracy,  resolution,  and 
convergence  of  the  numerical  solution.  On  the  windward  side 
of  the  plate  the  flow  is  essentially  a  stagnation-point-like 
flow  which  has  an  exact  solution  for  incompressible  flow. 
Figure  2  depicts  the  pressure  coefficient  at  the  windward  stag 
nation  point  as  a  function  of  the  initial  mesh  spacing  normal 
to  the  plate.  The  theoretical  pressure  coefficient  based  on 
isentropic  compression  at  Mo  =  0.5  is  Cp  =  1.0641.  Most  of 
the  results  for  this  figure  were  obtained  for  a  flow  domain 
substantially  smaller  than  the  final  domain.  The  outflow 
boundary  was  at  eight  semispans  downstream  of  the  plate  for 
the  first  set  of  four  runs,  labelled  points  (a),  (b) ,  (c) , 
and  (d) .  The  run  for  point  (e)  was  with  the  outflow  at 
twenty  semispans  downstream. 

The  first  case,  point  (a),  has  a  mesh  spacing  normal  to 
the  plate  of  0.001  (all  lengths  are  nondimensionalized  with 
respect  to  the  plate  semispan) .  The  mesh  point  distribution 
is  24  points  on  the  windward  side  of  the  plate,  21  points  on 
the  leeside,  and  78  points  parallel  to  the  plate  of  which  41 
points  are  on  the  plate  equally  spaced  at  Ax  =  0.025.  The 


stagnation  pressure  is  0.62  which  is  about  40%  lower  than  the 
isentropic  pressure  of  1.0641.  For  this  case  there  are  8 
points  within  the  boundary  layer  on  the  windward  side.  Thus 
there  are  only  16  points  between  the  inflow  boundary  and  the 
viscous  layer. 

Increasing  the  number  of  points  in  the  windward  inviscid 
region  at  the  expense  of  leeward  side  resolution  gives  the 
results  of  point  (b) .  The  stagnation  pressure  is  up  to  1.13 
which  compares  favorably  with  the  exact  value.  Further 
increases  in  the  number  of  mesh  points  in  the  windward  boundary 
layer  (to  12)  and  the  inviscid  region  (to  28  from  24)  and  pro¬ 
perly  resolving  the  leeward  region  gives  the  results  of  point 
(c) .  There  is  no  further  improvement. 

Reducing  the  mesh  spacing  normal  to  the  plate  to  0.0002 
(from  0.001)  and  increasing  the  number  of  boundary  layer  points 
to  15  gives  the  results  shown  by  point  (d) .  Here  the  pressure 
coefficient  is  1.67,  some  60%  too  high.  With  this  small  initial 
step  size,  the  mesh  must  be  rapidly  expanded  to  reach  the 
inviscid  layer  with  approximately  the  same  number  of  mesh  points 
The  thin  layer  approximation  used  in  the  code  is  not  correct 
for  highly  stretched  meshes  according  to  reference  4 .  The 
numerical  results  seem  to  verify  the  contentions  of  reference  4 . 
Since  it  was  not  possible  to  change  the  thin  layer  approxima¬ 
tion,  it  was  decided  to  choose  a  stretching  function  and  mesh 
spacing  so  that  the  thin  layer  approximation  as  used  in  the 
code  did  not  introduce  an  appreciable  error.  Thus  30  points 
were  selected  upwind  of  the  plate  with  8  points  in  the  boundary 
layer  and  an  initial  mesh  spacing  of  0.0005.  The  flow  field 
on  the  leeside  of  the  plate  with  the  outflow  boundary  at  8 
semispans  downstream  of  the  plate  indicated  that  the  separation 
bubble  was  not  closed.  Experimental  data  indicated  that  the 
separation  bubble  should  extend  approximately  8-10  semispans 
downstream  of  the  plate.  For  that  reason  the  outflow  boundary 


was  extended  to  20  semispans  downstream.  With  this  outflow 
boundary  48  points  were  needed  to  properly  resolve  the  separ¬ 
ated  region  and  the  shear  layers.  This  mesh  arrangement  is 
the  final  one  used  for  the  numerical  results  and  is  shown  as 
point  (e)  on  Figure  2. 

The  grid  network  for  this  final  configuration  is  shown  in 
Figures  3  and  4.  The  former  shows  the  entire  grid  and  the 
latter  shows  the  detailed  mesh  near  the  plate. 

The  convergence  history  of  the  numerical  solution  process 
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is  shown  in  Figure  5.  The  L  residual  error  drops  more  than 
3  orders  of  magnitudes  in  about  4000  iterations.  The  time 
step  was  taken  to  be  at  At  =  0.01  for  the  first  4000  iterations 
and  then  reduced  to  0.001  for  the  last  2000  iterations.  Since 
the  algorithm  is  implicit,  excessive  dissipation  occurs  for 
large  time  steps.  This  has  the  tendency  to  smear  the  shear 
layers  out  too  much.  By  reducing  the  time  steps  at  the  end  of 
the  run  the  shear  layers  can  be  "sharpened  up"  to  yield  more 
meaningful  results. 

The  final  check  to  verify  the  numerical  results  was  made 
by  comparing  the  boundary  layer  profiles  on  the  windward  side 
of  the  plate  with  the  exact  solution  available  for  the  incom¬ 
pressible  Navier-Stokes  equation  for  laminar  flows  (Rosenhead, 
ref.  5)  . 

A  brief  description  of  the  exact  solution  for  a  stagnation 
point  flow  is  appropriate.  A  stagnation  point  flow  on  a  flat 
plate  is  defined  as  a  flow  where  the  velocities  external  to 
the  boundary  are  given  by 

u„  =  kx  and  v_  =  -ky 
e  e 

where  k  is  a  constant  dependent  on  the  normalization  used  to 

nondimensionalize  u  .  v  ,  and  y  and  the  inviscid  flow  field. 
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The  subscript  'e'  refers  to  the  flow  external  of  the  boundary 
layer.  For  an  infinite  flat  plate  is  independent  of  y. 
However,  for  a  finite  flat  plate  ue  attains  a  maximum  at  the 
edge  of  the  boundary  layer  and  then  decays  with  increasing  y. 
This  is  shown  in  the  sketch  in  Figure  6.  Therefore,  for  the 


finite  flat  plate  u _  must  be  linear  in  x  for  the  flow  to 

max 


be  a  stagnation  point  flow.  This  linearity  is  shown  in 
Figure  6  for  the  numerical  solution.  As  shown  the  constant  k  is 


k  =  0.57  V/l 


The  exact  solution  (Rosenhead,  ref.  5)  is  given  in  terms  of  a 
similarity  solution,  and  the  similarity  variables  are 
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The  boundary  layer  profiles  in  terms  of  the  similarity 
variables  are  shown  in  the  next  two  figures.  The  numerical 
results  agree  fairly  well  with  the  exact  solution,  although  it 
is  apparent  that  the  boundary  layer  should  be  resolved  better 
with  more  mesh  points.  This  can  be  seen  by  the  difference 
between  the  numerical  and  exact  solution  in  Figures  7  and  8. 
However,  as  indicated  previously,  increasing  the  resolution 
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in  the  boundary  layer  at  the  expense  of  the  inviscid  region 
causes  errors  in  the  stagnation  pressure  and  reducing  the 
normal  step  size  near  the  plate  introduces  errors  from  the 
thin  layer  approximation.  So  the  mesh  point  distribution  as 
shown  in  Figures  3  and  4  was  accepted  as  a  workable  compro¬ 
mise.  However,  as  will  be  shown  later  the  thicker  numerical 
boundary  layer  on  the  windward  side  will  effect  the  sharpness 
of  the  shear  layer  coming  off  the  plate  edge. 

The  velocity  vector  plot  for  the  entire  flow  domain  is 
shown  in  the  several  parts  of  Figure  9.  The  details  near  the 
plate  are  in  Figure  10.  As  shown  there  is  a  large  separated 
flow  region  behind  the  plate  extending  some  11.5  semispans 
downstream.  Within  this  region  is  a  small  secondary  separated 
flow  region  just  inboard  of  the  leeside  of  the  plate  edge. 

This  can  be  seen  better  in  Figures  11  and  12.  The  boundary 
layer  profiles  are  also  readily  apparent  in  these  two  figures. 
Apparent  in  the  former  figure  is  the  constant  thickness  of  the 
boundary  layer  on  the  windward  side  of  the  plate .  Also  appar¬ 
ent  in  the  latter  figure  is  the  shear  layer  coming  off  the 
plate  edges. 

The  shear  layers  and  separated  flow  region  can  be  better 
observed  in  the  stream  function  contour  plots  shown  in  figures 
13,  14,  and  15.  These  three  figures  depict  the  stream  functions 
of  the  entire  flow  domain,  the  primary  separated  zone,  and  the 
secondary  separated  region,  respectively.  As  can  be  seen  the 
primary  separated  zone  extends  about  11.5  semispans  downstream 
of  the  plate.  The  secondary  bubble  although  small  does  extend 
over  20%  of  the  semispan  length  and  being  near  the  tip  of  the 
plate  has  an  influence  on  the  shedding  rate.  Contour  levels 
of  Mach  number,  pressure,  and  density  are  shown  in  Figures  16, 
17,  and  18. 
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2.4  Comparison  with  Experimental  Data 

An  extensive  set  of  experiments  was  conducted  on  two- 
dimensional  strakes  as  reported  in  reference  1.  In  this 
section  the  numerical  results  are  compared  with  those  experi¬ 
mental  data  and  some  theoretical  data. 

The  first  comparison,  shown  in  Figure  19,  is  the  static 
pressure  along  the  windward  side  stagnation  streamline.  The 
data  consists  of  measurements  on  the  plate  surface  plus,  in 
the  M  =  0.25  case,  static  pressures  deduced  from  velocity 
measurements  off  the  surface  using  the  incompressible  Bernoulli 
equation.  The  numerical  results  at  *  0.5  resemble  the 
behavior  of  the  experimental  data  for  =  0.25  and  would  be 
in  relatively  good  agreement  if  the  numerical  results  were 
displaced  downward  to  cause  the  Cp  on  the  plate  to  agree. 

The  static  pressure  on  the  face  of  the  plate  is  shown  in 
Figure  20.  Except  for  the  constant  shift  upward  of  the  numeri¬ 
cal  results  due  to  the  differences  of  pressures  at  the  stagna¬ 
tion  point,  the  numerical  and  experimental  data  agree  quite 
well.  We  note  that  the  experimental  pressure  coefficient  of 
1.06  at  the  stagnation  point  agrees  within  1%  of  the  theoreti¬ 
cal  isentropic  values  of  1.0641.  As  discussed  previously  the 
numerical  stagnation  pressure  depends  on  the  boundary  layer 
and  the  inviscid  region  resolution. 

There  are  two  principal  differences  between  the  numerical 
and  experimental  results.  One  is  the  sharpness  of  the  shear 
layer.  The  lack  of  mesh  resolution  has  the  effect  of  spreading 
the  numerical  shear  layer  and  thus  making  the  velocity  gradi¬ 
ents  less  severe  than  is  indicated  by  the  data.  Also  the  fact 
that  the  physical  model  has  a  cusp  while  the  numerical  model 
has  zero  plate  thickness  may  have  some  influence  on  the  differ¬ 
ences  in  the  gradients.  The  second  difference  is  due  to  the 


small  secondary  separation  bubble  near  the  leeward  edge  in  the 
numerical  solution  which  causes  the  v  component  to  be  every¬ 
where  downstream  (positive) ,  whereas  the  data  show  an  upstream 
flow  (negative  v)  in  the  recirculation  region  inside  the  shear 
layer.  The  effect  of  Reynolds  number  (0,15x10®  for  experiment 
and  0.5xio®  for  numerics)  is  not  important.  The  rapid  change 
of  the  axial  gradients  of  the  spanwise  velocity  at  the  same 
spanwise  location  is  shown  in  Figure  23. 

The  spanwise  variation  of  the  velocity  components  at  a 
station  just  behind  the  model  face  (y  =  0.2)  is  given  in 
Figures  24  and  25.  The  same  trend  is  shown  here,  i.e.  sharper 
experimental  shear  surfaces.  The  effect  of  the  secondary 
separation  bubble  is  apparent  in  the  numerical  results  for 
x  ^  1.  The  velocity  gradient  variations  at  the  same  axial 
station  are  shown  in  Figures  26  and  27.  The  vorticity  defined 
as 

3u  Sv 
5  ~  3y  *  3^ 

is  shown  in  Figure  28  as  a  function  of  span  at  the  same  axial 
location.  Figure  29  gives  corresponding  vorticity  transport 
rate  (vc)  variation.  Although  all  of  the  experimental  gradi¬ 
ents  are  higher  than  the  numerical  ones,  the  combination  of 
the  lower  gradients  and  higher  v  produce  vorticity  transport 
rates  similar  to  the  data. 

The  next  set  of  six  figures  (30  through  35)  presents  the 
same  sequence  of  results  as  shown  by  the  previous  six  figures 
except  at  an  axial  station  ^1.1  semispans  downstream  from  the 
model  face.  The  rapid  decay  of  the  spanwise  velocity  compo¬ 
nent  (u)  is  evident.  This  decay  is  much  more  rapid  for  the 
numerical  results.  This  is  due  not  to  any  physics  but  to  a 
lack  of  adequate  mesh  resolution.  At  this  location  of  the 
shear  layer  (y  -  1.1  and  x  =  1.7)  the  mesh  is  substantially 


coarser  than  at  the  plate  tip  (see  Figures  3  and  4) .  Since 

the  numerical  error  and  numerical  dissipation  is  proportional 
2  2 

to  0(Ax  ,  Ay  ,  AxAy)  for  the  implicit  scheme  used  in  the  code, 
the  increasing  coarseness  rapidly  increases  the  effective 
dissipation  of  the  scheme  and  thus  the  shear  layers  will  be 
spread  much  faster  than  dictated  by  the  physical  dissipation. 

The  lack  of  adequate  resolution  will  affect  the  other 
properties  of  the  shear  layer  as  shown  by  Figures  36,  37,  and 
38.  The  first  figure  is  the  axial  variation  of  the  stream- 
wise  vorticity  transport  rate,  R,  defined  by 
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R  =  J  v  C  dx 
o 

Shown  are  two  sets  of  experimental  data.  One  set  is  based  on 
the  above  formula  and  the  other  set  is  based  on  a  thin  shear 
layer  approximation  where 

3u  3v 
g  ~  3y  “  3x 

~  dv 
dx 

R  .  .  J  V  ||  dX  .  -  |  «<£> 
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and  o  and  i  represent  the  outer  and  inner  edges  of  the  shear 
layer,  respectively.  The  vorticity  transport  rate  was  computed 
from  the  exact  formula  by  the  trapezoidal  rule  for  the  numeri¬ 
cal  results.  As  can  be  seen  from  Figure  36,  the  numerical 
vorticity  transport  rate  is  low  compared  to  the  experimental 
results . 

Also  shown  on  this  figure  are  the  numerical  vorticity 
transport  rates  at  a  point  where  the  shear  layer  is  j  'st  coming 


off  the  plate  edge.  In  this  case  the  transport  rate  is 

6w 

R  =  -  |  u  C  dy 
-61 


where  u^w  is  the  spanwise  velocity  component  at  the  boundary 
layer  edge  on  the  windward  side  of  the  plate  tip  and  is  the 
corresponding  velocity  on  the  lee  side  of  the  plate  tip.  The 
calculated  vorticity  transport  rate  in  the  shear  layer  and  the 
boundary  layer  at  the  tip  agree  fairly  well.  The  differences 
between  the  numerical  calculation  and  the  experimental  data 
are  the  same  no  matter  how  they  are  computed. 

The  effect  of  decreasing  numerical  resolution  and  thus 
the  faster  spreading  rate  of  the  shear  layer  is  quite  evident 
in  this  figure  and  the  following  two  figures.  Figure  37 
depicts  the  integrated  vorticity  (y)  along  various  axial  sta¬ 
tions.  The  integrated  vorticity  is  defined  as 

00 

Y  =  j  ;  dx 
o 

which  becomes  with  the  thin  shear  layer  approximation 


Again  the  lower  values  for  the  numerical  results  are  evident. 
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The  circulation,  r,  can  be  computed  around  rectangular 
contours  at  many  y  locations,  each  contour  extending  spanwise 
from  the  symmetry  plane  to  the  far  field  boundary  and  having 
a  streamwise  length  of  10%  of  the  plate  semispan.  This  com¬ 
parison  is  shown  in  Figure  38.  The  circulation  is  essentially 
zero  upstream  of  the  plate  face,  a  rapid  increase  occurs  just 
after  the  face  to  a  final  peak  at  about  1.5  semispans  down¬ 
stream  of  the  face  and  then  a  gradual  decay.  As  expected  the 
numerical  results  are  lower  than  the  experimental  data. 

Another  difference  between  the  numerical  results  and  the 
experiment  is  the  lateral  extent  of  the  primary  separation 
bubble.  This  can  be  observed  by  comparing  figures  39  and  40 
with  Figure  23(b)  of  reference  1.  Shown  in  these  figures  are 
the  axial  velocity  profiles  at  various  axial  stations  behind 
the  model  face.  If  v  =  0  is  used  as  a  measure  of  the  size  of 
the  primary  separated  flow  region,  then  the  experiment  indi¬ 
cated  for  y  =  2.0  that  v  *  0  at  x  *  1.8  whereas  for  the  compu¬ 
tations  x  ~  1.0.  The  proper  measure  of  the  size  of  the  sepa¬ 
rated  region  is  the  dividing  streamline  (the  streamline 
emanating  from  the  tip  of  the  plate)  as  shown  in  Figure  15 
for  the  numerical  computations.  Unfortunately  this  is  not 
available  for  the  experiment.  It  can  be  approximated,  however, 
from  Figure  23(b)  of  reference  1  to  be  at  x  =  2.1  for  y  =  2.0. 
The  numerical  dividing  streamline  is  at  x  =  1.6  (Fig.  40). 

Thus  the  spanwise  extent  of  the  primary  separated  region  is 
much  smaller  for  the  numerical  results. 

There  are  at  least  two  possible  reasons  for  this  discre¬ 
pancy.  One  is  the  model  configuration.  The  experimental  model 
has  a  finite  thickness  plate,  finite  tip  angle,  and  a  splitter 
plate  whereas  the  numerical  configuration  consists  of  a  thin 
flat  plate  with  a  plane  of  symmetry  having  no  boundary  layer 
instead  of  the  splitter  plate  with  its  own  boundary  layer. 

These  differences  in  the  model  geometry  can  be  expected  to 


enlarge  the  separated  region.  The  second  reason  is  the  differ¬ 
ence  in  the  secondary  separation  regions.  In  the  numerical 
simulation,  there  is  a  secondary  separation  on  the  downstream 
side  of  the  plate  near  its  edge,  as  shown  on  the  left  sketch. 

The  experiment  has  a  secondary  separation  region 

on  the  splitter  plate  as  shown  on  the  right.  There  was  also  a 
difference  in  channel  half-width  (9  for  the  computation  vs.  12 
for  the  experiment)  although  this  should  have  very  little  effect 
on  the  problem. 

While  the  numerical  results  show  good  qualitative  agree¬ 
ment  with  the  experiment,  the  mesh  points  must  be  better  distri¬ 
buted  or  the  number  must  be  increased  to  properly  resolve  the 
shear  layers.  If  the  shear  layers  are  not  sufficiently  resolved 
then  little  accuracy  can  be  expected  from  the  nuemrical  compu¬ 
tations  about  the  vorticity  shedding  (transport)  rates  or  the 
circulation.  While  the  general  features  of  the  flow  field  are 
fairly  easy  and  inexpensive  to  compute,  accurate  details  of  the 
field  are  much  more  expensive  unless  flow  adaptive  meshes  are 
incorporated  into  the  numerical  computations. 

3.0  THREE-DIMENSIONAL  PROBLEM 

A  calculation  of  the  flow  field  about  a  75°/63°  double-delta 
wing  was  performed  on  the  CDC  7600  computer  using  the  version 
AIR3DP5  of  the  code  AIR3D  developed  by  Pulliam  and  Steger  and 
described  in  reference  6.  The  code  solves  the  three-dimensional 
"thin  layer"  Navier-Stokes  equations  by  an  implicit  approximate 
factorization  finite-difference  algorithm.  The  body  surface  must 
be  topologically  similar  to  a  hemisphere  cylinder. 


As  in  the  two-dimensional  code,  the  thin  layer  approxima¬ 
tion  means  that  the  only  second  derivatives  which  are  retained  in 
the  viscous  terms  are  those  normal  to  the  body.  This  is  necessi¬ 
tated  anyway  because  the  available  computer  storage  only  allows 
us  to  cluster  the  grid  in  one  direction,  namely  a  "boundary-layer 
like"  coordinate  5  normal  to  the  body. 

The  wing  configuration  used  in  the  calculations  was  that 
selected  originally  for  the  three-dimensional  wind  tunnel  test. 

The  sweep  of  the  forward  part  of  the  wing  (75°)  is  essentially 
that  of  an  aspect  ratio  1  delta  wing,  on  which  Hummel  (refs.  7 
and  8)  has  taken  surface  pressure  and  flow  visualization  data 
which  could  be  used  to  evaluate  the  computed  results  forward  of 
the  break  in  sweep.  The  sweep  of  the  aft  position  of  the  wing 
was  selected  first  to  be  not  so  different  than  that  of  the  for¬ 
ward  part  that  the  leading  edge  vortex  from  the  forward  part 
would  "tear  off"  at  the  break  in  sweep  and  second  to  be  close  to 
a  double-delta  wing  for  which  some  data  exist.  Wentz  (refs.  9  and 
10)  has  taken  data  on  a  number  of  double-delta  wings,  one  of  which 
has  sweeps  of  75°  and  62°.  The  data  consists  of  forces  and 
moments,  oil  flow,  surface  pressures,  and  velocities  above  the 
wing.  This  configuration  is  sufficiently  close,  to  the  750/63° 
wing  to  permit  evaluation  of  the  numerical  results. 

3.1  Body  Geometry  and  Mesh  Arrangement 

The  computational  body  is  pictured  in  orthographic  projec¬ 
tion  in  Figure  41.  The  planform  is  a  double-delta  shape  with 
sweep  angles  of  75°  and  63°,  with  the  break  in  sweep  occurring  at 
midchord  (x=0.5).  Because  of  coding  requirements,  the  body  is 
extended  aft  of  the  root  chord  position  at  constant  semispan  (the 
"sting") ,  thus  implying  a  second  break  in  sweep  to  90°  at  the  tail 
plane  x  =  1.0.  All  dimensions  are  normalized  to  root  chord. 

The  cross  sections  are  ellipses  chosen  so  that  the  vertical 
axis  is  .01  root  chord  and  the  horizontal  axis  is  equal  to  the 
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semispan.  This  sectional  shape  allows  the  use  of  orthogonal 
elliptic  coordinates  in  the  crossflow  planes,  as  can  be  seen  in 
the  grid  plots  (Fig.  42) .  Since  the  vertical  axis  is  constant  and 
the  semispan  varies,  the  axis  ratio  e  of  the  ellipses  is  not  con¬ 
stant  but  varies  from  more  than  one  near  the  nose  to  about  .013 
at  and  aft  of  x  =  1 .  Hence,  the  configuration  is  not  conical. 

In  order  to  avoid  a  coordinate  singularity  at  the  nose, 
the  body  was  truncated  at  x=  .01  and  a  prolate  spheroidal  "cap" 
was  appended  thereto. 


Physical  coordinates  (x,y,z)  are  chosen  so  that  x  increases 
from  zero  along  the  chord  line  going  aft,  y  increase  toward  star¬ 
board  and  z  increases  upward.  Computational  coordinates  (£,n»£) 
are  chosen  so  that  z,  =  0  is  the  body  surface.  Except  on  the  cap, 
£  measures  the  axial  distance  and  is  equal  to  x.  n  is  an  azi¬ 
muthal  angle  and  z,  increases  radially  outward.  Both  systems  of 
coordinates  are  right-handed.  On  the  cap,  coordinate  £  measures 
the  polar  angle. 


Let  b(£),  £  >  .01,  be  the  semispan  as  a  function  of  radius. 

Thus , 


b  ( £) 


£  tan  15°,  .01  <  £  <  0.5 

0.5  tan  15°  +  (£  -  0.5)  tan  27°,  0.5  <  £  <  1 
0.5  (tan  15°  +  tan  21°),  1.0  <  £ 


Let  b  .  =  .01  tan  15°  and  b  =  .5  (tan  15°  +  tan  21°).  Then, 

min  max 

the  transformation  between  physical  and  computational  coordinates 


is  as  follows: 


On  the  cap: 

x  =  .01  -  /  c2  +  b2 .  cos  n  cos  £ 

min 

y  =  /  £2  +  b2.  cos  n  sin  £ 
min 

z  -  z,  sin  n  sin  £ 


!» 


27 


On  the  wing: 


x  =  5 

y  =  /  ;2  +  It>  ( C )  3 2  cos  n 
z  =  ;  sin  n 


On  the  sting: 


x  =  5 

y  =  /c2  +  b2  cos  n 
max 

z  =  £  sin  n 

The  computation  takes  place  in  a  cubical  domain  pictured 
in  Figure  43  (see  also  Fig.  2(b)  of  ref.  1).  The  computational 
coordinates  are  discretized  so  that 


K  =  C(J) /  J  =  1/30 
n  *  n (l) ,  l  =  1,30 
Z,  ■  C(K) ,  K  *  1,24 


The  cross  flow  planes  are  represented  by  J  *  constant  computa¬ 
tional  planes,  and  are  tabulated  in  Table  1.  The  nose  is  at 
J  =  8  and  the  sting  begins  at  J  =  24.  There  is  a  singular 
line  at  J  -  K  «  1  corresponding  to  the  pole  of  the  cap.  This 
singularity  is  avoided  in  the  finite-difference  version  of 
the  problem  because  the  fluxes  in  the  £-  and  C-  directions 
vanish  identically  on  this  line  (ref.  6)  . 

The  mesh  is  clustered  in  the  axial  direction  near  the 
nose,  break  in  sweep,  and  root  chord,  and  in  the  circumferential 


TABLE  I.  75*763°  DOUBLE  DELTA  WING  CALCULATION 


Mesh  Crossflow  Planes 


Index  Chord 
J  z(J) 


Semispan 
b  (J) 


Ellipse 
Axis  Ratio 


Nose  region  "cap' 
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Indicat 


.0100 
.0341 
.0824 
.1579 
.2496 
.  3408 
.4162 
.4683 
.50 
.  5025 
.  5384 
.5936 
.6712 
.7632 
.8543 
.9287 
.9761 
1.00 
1.031 
1.123 
1.389 
2.166 
4.425 

11.000 

es  plots 


.0026 

.0091 

.0220 

.0423 

.0669 

.0913 

.1115 

.1255 


■Break 


.1353 
.  1535 
.  1816 
.2212 
.  2681 
.3145 
.3524 
.3765 
.38874 


1.92 


.118 

.075 

.055 

.045 


in  sweep 
.037 
.0326 
.0275 
.0226 
.0187 
.0159 
.0142 
.0133 
.0129 


2.547 

2.485 

2.360 

2.166 

1.930 

1.696 

1.502 

1.368 


1.280 
1 .187 
1.045 
.846 
.609 
.375 
.183 


-  .080 


-1.001 
-2.999 
-8.811 
25. 724 


w 


direction  near  about  80%  semispan  on  the  leeward  side  where  we 
expected  that  the  primary  leading  edge  separation  vortex  would 
form.  In  the  circumferential  direction,  one  extra  point 
beyond  the  plane  of  symmetry  at  each  boundary  is  required  by 
the  code  for  implementing  the  boundary  conditions.  The  radial 
direction  is  clustered  so  that  8  of  the  24  points  lie  within 
.15  chordlengths  of  the  body  surface.  The  grid  is  exponen¬ 
tially  stretched  out  to  an  outer  boundary  10  chord-lengths 
away.  Figures  42  and  44  give  typical  views  of  the  grid. 

The  grid  generation  routine  calculates  x(J,K,L) ,  y(J,K,L), 
z ( J ,K ,L)  for  each  of  the  21,600  grid  points  and  stores  these 
on  tape  for  use  in  the  calculation  itself  and  the  graphics 
display  package. 

3.2  Calculation  Procedure 

The  free-stream  Mach  number  of  the  calculation  was  0.5. 
This  was  large  enough  to  let  the  calculation  converge  but  small 
enough  to  prevent  the  development  of  a  large  supersonic  region 
near  the  edges.  The  angle  of  attack  was  15°,  so  the  Mach 
number  in  the  crossflow  plane  was  0.13.  The  Reynolds  number 
was  0.5*106  based  on  root  chord,  and  a  laminar  boundary  layer 
was  selected.  The  Prandtl  number  was  0.72  and  the  free-stream 
temperature  was  500 °K. 

The  calculations  were  run  on  a  CDC  7600  computer  at  Ames 
Research  Center.  Two  smoothing  parameters  were  used  in  AIR3D, 
as  described  in  reference  6.  The  explicit  smoothing  parameter 
was  set  at  0.1  and  the  implicit  parameter  was  set  at  1.0, 
based  on  previous  experience  with  the  code . 

Following  past  experience  with  the  code,  the  boundary 

conditions  were  entered  during  the  first  30  iterations,  and 

-4 

the  code  was  run  at  a  timestep  of  2 ><10  ,  corresponding  to 

a  Courant  number  of  about  4  for  500  iterations.  The  timestep 
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was  increased  to  5x10  (CFL  =  10)  after  100  iterations  and 
-4 

to  10x10  (CFL  :  20)  after  1400  iterations.  The  program  was 
run  for  9500  iterations,  until  the  residual  was  down  by  about 
a  factor  of  30.  The  CPU  time  per  iteration  was  about  17 
seconds.  Thus,  even  on  this  relatively  coarse  grid,  the  calcu¬ 
lation  required  more  than  80  hours  of  CPU  time.  The  convergence 
history  is  shown  in  Figure  45. 


3.3  Results 

The  flow  field  (i.e.,  values  of  the  density,  three  compo¬ 
nents  of  the  velocity,  and  the  internal  energy  at  every  grid 
point)  was  written  into  a  disc  file  for  analysis  and  graphic 
display. 

The  closest  experimental  comparison  to  this  calculation 
is  the  work  of  Wentz  and  McMahon  (refs.  9  and  10)  who  used  a 
75°/62°  double  delta  pictured  in  Figure  46.  There  are  two 
primary  differences  between  the  experiment  and  the  calcula¬ 
tion.  The  experimental  model  has  a  fuselage  attached,  and 
has  a  sharp  trailing  edge,  while  the  computational  model  has 
neither.  Also,  the  experiment  was  run  at  M^  ~  .16  while  the 
calculation  was  at  M  =  .50  (the  latter  was  necessitated  in 

oo 

order  to  obtain  convergence  within  a  reasonable  time).  Wentz's 
most  complete  data  were  obtained  at  angles  of  attack  of  10°  or 
20°,  so  the  comparisons  at  15°  are  somewhat  limited. 

Plots  of  the  flow  field  were  obtained  for  every  crossflow 
plane  aft  of  the  nose.  Four  sets  of  these  are  reproduced  here. 
Figures  47 (a) -(h)  show  the  flow  field  at  about  one-quarter 
root  chord,  halfway  between  the  nose  and  the  break  in  sweep. 
Figures  48 (a) -(h)  show  the  flow  field  just  aft  of  the  break  in 
sweep  at  mid-chord.  Figures  49 (a) -(h)  show  the  flow  field  at 
three-quarters  chord  and  Figures  50 (a) -(h)  show  the  flow  field 
at  the  root  chord  where  the  sting  begins. 
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For  each  of  these  chordwise  stations  are  plotted  (a)  the 
mesh,  (b)  an  overall  view  of  the  crossflow  velocity  field, 

(c)  a  closeup  of  the  velocity  field  near  the  tip,  (d)  a  plot 
of  Cp  along  the  wing  surface,  and  (e)  contours  of  the  cross- 
flow  Mach  number,  with  20  equally  spaced  levels  between  maxi¬ 
mum  and  minimum  Mach  number,  AM  *  .05.  Figure  (f)  of  each 
set  shows  the  vorticity  contours,  with  levels  exponentially 
concentrated  near  the  maximum  value.  In  all  cases  the  con¬ 
centration  is  very  near  the  tip.  Figure  (g)  of  each  set 
shows  contours  of  density  and  Figure  (h)  of  each  set  shows 
contours  of  total  pressure,  with  equal  increments  between 
contours . 

The  most  important  qualitative  features  of  the  flow  can 
be  seen  in  the  velocity  vector  plots.  One  can  see  the  primary 
vortex  system  forming  at  about  x  =  .25.  The  center  of  this 
vorticity  moves  inboard  until  at  the  tail  plane,  its  center 
is  at  about  y  =  .22  and  z  =  .12,  or,  in  terms  of  semispan, 
y/b  «  .6  and  z/b  s  .3.  This  may  be  compared  with  the  simpli¬ 
fied  theory  of  Brown  and  Michael  (ref.  4)  which  predicts  a 
location  y/b  =  .88,  z/b  =  .14.  We  can  see  that  the  onset  of 
primary  separation  occurs  at  about  x  =  .25  and  that  the  separ¬ 
ation  line  reaches  the  plane  of  symmetry  at  about  x  =  .5. 

One  can  see  the  development  of  secondary  separation  near  the 
tip  in  Figure  50(c),  but  the  mesh  is  too  coarse  to  resolve  it. 
In  the  same  figure,  note  the  typical  boundary  layer  profile  on 
the  lower  surface.  This  feature  is  qualitatively  correct,  but 
the  profile  should  have  a  much  smaller  thickness.  This  is 
again  due  to  the  inadequacy  of  the  grid,  even  though  it  is  the 
largest  that  can  be  run  on  the  CDC  7600. 

Looking  at  the  pressure  contours,  we  find  that  the  Cp 
distribution  is  very  sharply  peaked  near  the  tip,  unlike  the 
experimental  (i.e.,  ref.  8  for  a  delta  wing  and  ref.  9  for  a 
double  delta)  results  which  show  milder  peaks  at  about  70% 
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semi span.  Indeed,  there  is  a  small  region  where  the  crossflow 
is  supersonic,  even  though  the  crossflow  freestream  Mach 
number  is  about  0.13.  Perhaps  this  is  due  in  part  to  the 
unevenness  of  the  mesh  in  the  region  and  a  large  concentra¬ 
tion  of  grid  points  near  the  tip. 

The  other  plots  are  more  or  less  self-explanatory.  Con¬ 
tours  of  the  pressure  coefficient  on  the  wing  surface  are  pre¬ 
sented  in  Figure  51.  Wentz's  experimental  result  is  presented 
in  Figure  52.  Several  features  are  different.  In  the  calcu¬ 
lation,  we  find  a  minimum  Cp  of  -2.45,  which  corresponds  to  a 
slightly  supersonic  Mach  number.  Wentz's  minimum  Cp  is  about 
-1.95.  Part  of  this  is  due  to  compressibility.  Applying  the 
Prandtl-Glauert  correction  (the  effect  of  compressibility  to 
first  order  -  see  ref.  11)  yields  the  following  comparison: 

Wentz's  value  (at  M  =  .16)  is  corrected  to 
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1.93 


and  our  minimum  (at  M  =  .5)  is  corrected  to 
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2.12 


The  predicted  vortex  lift  is  more  concentrated  and  further 
outboard  (near  95%  semispan)  than  the  experimental  result, 
which  shows  the  locus  of  minimum  Cp  at  roughly  70%  of  the  semi¬ 
span. 

On  the  computational  wing  there  is  not  as  much  lift  on 
the  portion  near  the  nose.  This  is  due,  at  least  in  part,  to 
the  absence  of  a  forebody  which  for  the  Wentz  results  causes 
additional  upwash  near  the  wing  leading  edge  where  the  span 
is  the  same  order  as  the  body  radius. 


33 


The  total  circulation,  nondimensionalized  as  F  *  T/2V  b 

00  max 

to  compare  with  Wentz's  measurements,  was  calculated  over  half 
the  wing  by  performing  the  contour  integration  |  $*d t  over  the 
contours  shown  in  Figure  53.  Contour  encloses  the  entire 
flow  field  bounded  by  the  plane  of  symmetry.  Thus  the  contri¬ 
bution  on  the  wing  surfaces  vanishes  due  to  the  boundary 
condition  v  *  0.  Contour  Cv  encloses  the  upper  quadrant  of 
the  flow  field  but  excludes  the  first  three  grid  points  from 
the  wing  surface  (i.e.  the  "boundary  layer").  This  is  a  closer 
approximation  to  the  contour  Wentz  used,  since  he  could  not 
obtain  measurements  in  the  boundary  layer.  The  circulation 
is  presented  as  a  function  of  axial  distance  in  Figure  54. 
Theoretical  values  according  to  Jones  (linear  theory,  ref.  3) 
and  Brown  and  Michael  (concentrated  vortex  core,  ref.  12)  are 
shown  for  a  =  15°.  Wentz's  experimental  results  for  a  =  20° 
and  his  contour  are  shown  in  Figure  55.  In  both  cases  the 
results  lie  between  the  two  theoretical  curves.  The  falloff 
in  circulation  over  the  aft  half  of  the  wing  is  due  to  the 
presence  of  the  trailing  edge  in  the  experiment  and  the  "sting" 
in  the  calculation.  A  more  recent  experiment  by  Sforza  and 
Smorto  (refs.  14-15)  also  has  shown  a  value  for  the  circulation 
which  lies  between  the  Jones  and  Brown-Michael  theories. 

Integrating  the  pressure  coefficient  over  the  body  surface, 
we  find  a  normal  force  coefficient 


and  a  lift  coefficient  CT  =  C.,cos(15°)  =  0.747.  Here  S  is  the 

1  L  N 

planform  area  and  q  =  j  pU  is  the  dynamic  pressure.  Normalized 
by  aspect  ratio  we  obtain  (C^/AR)  =  0.420  for  an  aspect  ratio 
of  1.84.  Wentz  (refs.  9  and  10)  obtains  the  values  CL  =  0.70 
and  (C^/AR)  *  0.43,  respectively,  his  aspect  ratio  being  1.61. 
While  the  agreement  is  very  good,  it  is  difficult  to  interpret 


this  result,  in  that  the  nonlinear  portion  of  the  lift  can  be 
considered  the  product  of  the  vortex  strength  and  span,  which 
could  both  be  somewhat  in  error  while  still  producing  the 
correct  overall  lift. 

Again  by  integration  we  obtain  a  center  of  pressure  of 

x  *0.77  compared  to  Wentz's  x  =0.58.  That  is  to  say, 
cp  cp 

the  calculation  shows  rather  more  lift  on  the  aft  sections. 

This  is  most  likely  attributable  to  the  presence  of  the  sting 
on  the  wing,  which  maintains  the  lift  near  the  "trailing  edge", 

as  can  be  seen  from  the  C  contour  plots . 

P 

3.4  Notes  on  the  ILLIAC  IV  Calculation 

The  plan  in  this  work  was  to  run  the  AIR3D  code  on  the 
CDC  7600  to  evaluate  its  ability  to  treat  slender  thin  wings 
and  if  successful  then  move  to  the  ILLIAC  IV  where  a  finer 
mesh  could  be  used.  The  calculation  performed  on  the  CDC  7600, 
while  generally  successful,  did  not  permit  a  sufficient  number 
of  grid  points  to  capture  the  boundary  layer.  A  typical 
boundary  layer  velocity  profile  occurs  on  the  windward  side 
[Fig.  50(c)],  but  it  is  "stretched  out"  to  about  ten  times 
what  one  would  normally  expect  for  a  Reynolds  number  of  500,000. 
Clearly,  it  was  necessary  to  calculate  on  a  finer  mesh.  Thus, 
the  calculation  was  programmed  for  the  ILLIAC  IV  computer. 

Fortunately,  the  program  AIR3D  has  been  implemented  on  the 
ILLIAC  IV  and  run  for  simple  geometries  (ref.  17).  The  capacity 
of  the  ILLIAC  IV  in  64-bit  mode  is  80  axial  *  56  circumferential 
x  48  radial  grid  points,  a  total  of  215,040  or  about  ten  times 
what  could  be  done  on  the  CDC  7600.  Simpler  geometries,  like 
the  hemisphere  cylinder,  has  been  run  before,  but  this  new 
calculation  would  be  more  complicated,  and  on  a  slightly  larger 
mesh  than  had  been  run  before. 


The  unique  nature  of  the  ILLIAC  hardware  requires  equally 
unique  data  structures.  Thus,  there  was  a  time-consuming  pro¬ 
cess  of  generating  the  mesh,  and  formatting  it  and  the  data 
onto  files  which  could  be  used  by  the  ILLIAC. 

A  new  problem  was  chosen,  with  a  somewhat  simpler  geometry 
in  order  to  concentrate  on  modeling  the  leading  edge  vortex  to 
see  if  better  agreement  could  be  obtained  on  spanwise  location 
and  surface  pressure  distribution.  The  problem  selected  was 
a  simple  delta  wing  of  unit  aspect  ratio,  at  20.5°  angle  of 
attack,  such  as  had  been  investigated  by  Hummel  (refs.  7  and  8) 

It  was  clear  from  experience  with  the  CDC  7600  calculation 
that  the  orthogonal  mesh  used  there  inhibited  convergence. 

Thus  a  new,  clustered  mesh,  nonorthogonal  in  the  crossflow 
planes,  was  generated  for  the  new  calculation. 

The  new  mesh  was  written  onto  tape  and  the  code  was 
assembled  on  tapes  for  the  ILLIAC  IV  computer.  The  code  was 
run  for  100  iterations,  long  enough  to  make  sure  it  was  running 
properly  and  to  impose  the  boundary  conditions.  The  CPU  time 
was  about  72  seconds  per  iteration.  Unfortunately,  the  ILLIAC 
IV  was  decomissioned  before  further  iterations  could  be  run, 
and  no  further  work  was  possible. 

The  code  AIR3D  is  being  rewritten  for  use  on  the  CRAY-1 
but  it  will  not  be  possible  to  run  the  code  on  the  CRAY  until 
the  size  of  the  memory  is  extended  to  8  million  words,  which 
will  not  occur  until  1983. 

3.5  Note  on  Mesh  Generation 

It  was  noted  during  the  work  that  the  largest  permissible 
timestep  for  a  given  Courant  number  was  controlled  by  the  mesh 
cell  size  near  the  leading  edge.  This  is  illustrated  in  the 
closeup  of  the  mesh.  Figure  56.  Consider  a  crossflow  plane 
where  the  semispan  is  b,  and  the  minimum  radius  at  the  edge 
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has  the  value  o.  Then  the  C  =  constant  surfaces,  those  surfaces 
parallel  to  the  wing  surface,  are  given  by  the  ellipses. 


=  Ky  + 


a)2  +  b2]1/2 


Differentiating,  we  find 


z 


+  b2 


For  small  z,  we  have  dy  :  g  (dz) .  Since  z/b  is  of  the  order 
0.01  near  the  z  =  1  plane,  this  compression  of  the  z,  =  const, 
surfaces  near  the  edge  slowed  down  the  maximum  timestep  by  at 
least  an  order  of  magnitude. 

A  remedy  for  this,  which  was  checked  in  working  out 
the  mesh  for  the  ILLIAC  IV  calculation  was  to  use  the  non- 
orthogonal  mesh 


z  =  [(y  +  o)n  +  bn]1/n 

which  is  orthogonal  for  n  =  2.  For  n  =  1,  there  were  problems 
with  skewness  of  the  mesh;  n  =  1.5  provided  a  good  compromise, 
and  allowed  a  tenfold  larger  time-step  for  a  given  Courant 
number.  The  visible  improvement  over  Figure  56  is  seen  in 
Figure  57,  where  the  same  region  of  he  edge  is  shown  in 
the  new  mesh. 


4.  CONCLUDING  REMARKS 

The  objective  of  this  portion  of  the  strake  vorticity 
shedding  work  was  to  demonstrate  that  computational  solutions 
could  be  used  to  enhance  one's  knowledge  of  an  experimentally 
determined  flow,  particularly  for  regions  in  the  flow  where 
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measurements  are  very  difficult  to  obtain.  In  order  to  use  a 
computational  solution  in  this  manner,  the  important  features 
of  the  computational  solution  must  agree  well  with  those  of 
the  experiment.  In  the  two  cases  considered  here,  the  flow 
is  dominated  by  shear  layers  separating  from  an  edge  which 
depends  strongly  on  the  boundary  layer  characteristics  of  the 
flow  approaching  the  edge.  Thus,  it  is  important  to  t  ■»  able 
to  adequately  resolve  and  describe  in  the  computation  the 
boundary  layers  on  the  surfaces  approaching  the  edge  and  the 
large  flow  property  gradients  in  the  shear  layer  adjacent  to 
the  edge . 

The  two  cases  considered  here  were  a  two-dimensional  flow 
normal  to  a  thin  flat  plate  and  a  three-dimensional  flow  over 
a  thin  slender  wing  at  moderate  angle  of  attack.  In  spite  of 
the  fact  that  the  two  cases  are  conceptually  simple,  there  is 
little  flow  field  data  available.  In  the  two-dimensional  case, 
the  data  consist  of  laser  velocimeter  measurements  taken  as  an 
earlier  part  of  this  study  and  data  taken  6  decades  ago  with  a 
hot  wire  anemometer.  In  the  three-dimensional  case,  there  are 
flow  field  data  above  the  wing  and  pressures  on  the  wing,  but 
no  data  near  the  edge  to  define  the  shear  layer  and  vorticity 
shedding  rate. 

It  is  important  to  note  that  when  the  present  investiga¬ 
tion  was  initiated,  neither  the  two-dimensional  nor  the  three- 
dimensional  computer  code  had  been  used  to  compute  flows  with 
substantial  separation.  Hence  the  present  work  can  be  viewed 
as  a  preliminary  study  of  the  ability  of  these  computer  codes 
to  capture  and  represent  the  various  separated  flow  phenomena. 

In  both  cases,  the  important  flow  features  of  the  compu¬ 
tational  flow  agreed  qualitatively  with  those  shown  by  the 
data.  Specifically,  the  computational  boundary  layers  on  the 
windward  and  leeward  surfaces  near  the  edge  generated  vorticity 
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which  appears  in  a  shear  layer  leaving  the  edge.  In  the  two- 
dimensional  case,  the  primary  separation  bubble  behind  the 
plate  plate  had  nearly  the  same  axial  and  lateral  extent  as 
shown  by  the  experiment.  In  the  three-dimensional  case,  the 
computational  vorticity  appeared  in  the  flow  above  the  wing, 
although  it  did  not  tend  to  roll  up  into  a  concentrated  vortex 
as  the  experiment  shows  [compare  Figs.  48(b)  and  58], 

In  both  cases,  computational  mesh  size  and  arrangement 
limitations  are  an  important  factor  in  the  comparisons  with 
data.  It  was  not  possible  to  keep  enough  mesh  points  in  the 
boundary  layer  and  still  define  well  enough  the  flow  away  from 
the  surface.  This  was  particularly  so  for  the  three-dimensional 
case,  which  had  a  much  more  severe  storage  problem  for  the 
computer.  This  is  felt  to  be  a  major  limitation  in  achieving 
the  objective  of  using  computational  flow  fields  to  augment 
experimental  measurements. 

In  both  cases,  the  computational  problem  was  slightly, 
but  probably  importantly,  different  than  the  experimental 
problem.  In  the  two-dimensional  case,  the  experimental  flat 
plate  normal  to  the  flow  was  required  to  have  a  splitter  plate 
downstream  to  enforce  steady  flow  and  a  plane  of  symmetry. 

Thus,  the  boundary  condition  on  the  plane  of  symmetry  in  the 
experiment  was  a  no-slip  zero-velocity  condition  yielding  a 
forward  moving  boundary  layer  which  caused  secondary  separation 
to  occur  on  the  splitter  plate  just  downstream  of  the  plate. 

The  computation  had  a  plane  of  symmetry  condition  only,  with 
a  no-slip  condition  only  on  the  back  face  of  the  plate,  which 
yielded  secondary  separation  on  the  back  face  of  the  plate  near 
the  edge.  These  differences  were  probably  responsible  for  some 
of  the  differences  in  the  vorticity  shedding  rate  and  location 
of  the  shear  layer. 
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In  the  three-dimensional  case,  the  forward  portion  of  the 
wing  was  not  thin  because  of  a  limitation  in  the  "body"  geometry 
the  code  could  handle.  Also  the  code  could  not  treat  a  trailing 
edge.  Both  of  these  factors  would  affect  the  load  distribution, 
as  compared  to  the  experiment  and  possibly  the  general  charac¬ 
ter  of  the  vorticity-dominated  region  above  the  wing. 

In  the  three-dimensional  case,  there  was  also  a  qualita¬ 
tive  difference  between  the  computational  and  experimental  flow. 
The  computational  suction  peak  was  further  outboard  than  the 
measured  one,  and  the  computational  vorticity  appears  to  be 
more  diffuse  above  the  wing  than  the  data  would  indicate.  The 
extent  to  which  these  differences  are  dependent  on  the  mesh 
size  and  arrangement  and  the  wing  geometry  model  is  not  clear. 

In  the  three-dimensional  calculation,  it  became  clear  that 
the  design  of  the  coordinate  mesh  has  a  significant  effect  on 
the  convergence  time  of  the  solution,  as  well  as  on  the  accu¬ 
racy  with  which  the  features  of  the  flow  can  be  predicted  with 
a  fixed  mesh  (in  which  the  mesh  points  do  not  move  as  the  solu¬ 
tion  progresses) ;  one  must  pre-cluster  the  points  where  the  flow 
features  are  expected  to  appear.  Here,  we  clustered  the  points 
near  the  vortex  loci  shown  by  the  data,  but  the  code  caused 
the  vorticity  to  occur  at  a  different  spanwise  location.  In 
an  adaptive  mesh  algorithm,  the  points  are  allowed  to  move  in 
such  a  manner  as  to  "capture"  the  flow  phenomena.  Although 
the  storage  problem  is  not  eased,  an  adaptive  mesh  algorithm 
should  provide  better  definition  of  the  flow,  and  thus  an 
improved  solution  for  a  given  amount  of  computer  time. 

In  summary,  it  appears  that  for  two  dimensional  flows  with 
separation,  there  is  probably  a  reasonable  expectation  for 
using  computational  solutions  for  augmenting  experimental  ones, 
provided  care  is  exercised  in  modeling  the  flow  and  generating 
the  mesh.  For  three-dimensional  flows,  the  problem  appears 
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more  difficult  because  of  the  mesh  requirements  and  the  addi¬ 
tional  computer  storage  due  to  the  third  dimension  and  would 
probably  require  the  development  of  an  adaptive  mesh  algorithm. 
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Figure  1.  Two-dimensional  flow  domain 
semispan  length  It  =  1 . 
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Figure  2.  Windward  stagnation  pressure  coefficient 
as  a  function  of  mesh  spacing  near  plate  DY  and 
number  of  ncints  in  bound a rv  Isvcr. 
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Figure  4.  Detail  of  grid  plot  near  the 
2-dimensional  strake. 
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External  velocity  of  stagnation  point  flows 
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Figure  7.  Stagnation  point  boundary  layer  velocity  profile 
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(a)  Entire  flow  domain 
Figure  9.  Velocity  vector  plot 


(c)  Flow  detail  just  downstream  of  plate 
Figure  9.  (Continued). 
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{ e )  Flow  detail  well  downstream  of  plate 
Figure  9.  (Concluded). 
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Figure  10.  Detail  of  velocity  vector  plot 
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(a)  Inner  part  of  plate 

Figure  11.  Detailed  velocity  vector  plot  of  the 
flat  plate  boundary  layer. 
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(b)  Region  near  plate  edge 
Figure  11.  (Concluded). 


Figure  12.  Detailed  velocity  vector  plot  of  secondary 
separation  bubble  and  shear  surface 
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Figure  15.  Stream  function  contours  of  secondary 
separation  bubble  and  shear  surface. 


Figure  17.  Pressure  contours  of 
plate  at  Mo  =  0.5  and  Re  =  0. 
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Figure  19.  Windward  stagnation  line  static  pressure 

coefficient  (x/£  =  0). 
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Figure  22.  Axial  variation  of  spanwise  velocity 

at  x  =  1.3. 
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Figure  27.  Axial  gradient  of  spanwi se  velocity  at  y 


{ 


5 


1  2  3  4  5  6 

X 

30.  Axial  velocity  profile  at  y  =  1.1. 
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Figure  36. 


Comparisons  of  vorticity  transport  rates 


Figure  41.  Wing  configuration 
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re  42.  Crossflow  plane  elliptic  coordinates. 
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iguration  tested  by  Wentz  and  McMahon 


(b)  Crossflow  plane  velocity  vectors 
Figure  47.  (Continued) . 
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(c)  Crossflow  plane  velocity  vectors 
Figure  47.  (Continued). 
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Figure  47.  (Continued). 
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75/63  DOUBLE  DELTA  J  =  12  Z=.2496  B=.0669 


(h)  Pressure  contours 
Figure  47.  (Concluded) 
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(d)  Plot  of  Cp 
Figure  48.  (Continued) 


(g)  Density  contours 
Figure  48.  (Continued) 
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(h)  Pressure  contours 
Figure  48.  (Concluded) 
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75/63  DOUBLE  DELTA  J=20  Z=.7632  B=.2681 


(b)  Crossflow  plane  velocity  vectors 
Figure  49.  Continued 
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(f)  Vorticity  contours 
Figure  49.  (Continued) 
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(g)  Density  contours 
Figure  49.  (Continued) 


75/63  DOUBLE  DELTA  J=20  Z=.7632  B=.268i 


M„-  0.S0  o-  9000  SO- 10* 


(h)  Pressure  contours 
Figure  49.  (Concluded). 


117 


75/63  DOUBLE  DELTA  J=24  Z=l.  B--.3887  TAIL  PLANE 
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(b)  Crossflow  plane  velocity  vectors. 
Figure  50.  (Continued) . 
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(c)  Crossflow  plane  velocity  vectors 
Figure  50.  (Continued). 
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(e)  Crossflow  Mach  number  contours 
Figure  50.  (Continued) . 
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75/63  DOUBLE  DELTA  J=24  Z=l.  B=.3887  TAIL  PLANE 


(g)  Density  contours 
Figure  50.  Continued. 
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75/63  DOUBLE  DELTA  J=24  Z=l.  B=.3887  TAIL  PLANE 


(h)  Pressure  contours 
Figure  50.  (Continued). 


Figure  55.  Comparison  of  measured  and  predicted  circulation 
results  from  Wentz.  (Ref.  10). 
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Figure  57.  Close  up  of  nonorthogonal  mesh  near  edge. 
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